*************************************************************
*** Monte Carlo experiments (SAR DGP) for the WWW JoP Project
***
*** Created: 9-19-14
*** Modified: 9-5-19
***
*************************************************************

clear
clear matrix
clear mata
set maxvar 32000

version 14

do "Programs\SAR Programs.do"

*************************************************************
*************************************************************
*** SAR Data-Generating Process
*************************************************************
*************************************************************


*************************************************************
*** Set up the parameters of the MC experiments
*************************************************************
scalar b1 = 1				/* Parameter for beta_x1 */
matrix I = I(200)			/* Identity matrix with N observations */

local run = 1				/* Run number (for running multiple do files at once) */
local s = 2					/* Standard deviation of the errors */
local MCsims = 1000			/* Number of simulations in Monte Carlo experiments */
local QIsims = 1000			/* Number of simulations for uncertainty in quantities of interest */
local obs = colsof(I)		/* Number of observations */
local seed = 648			/* Set the seed */

set seed `seed'
set obs `obs'
if `obs' > `QIsims' {
	set matsize `obs'
}
else {
	set matsize `QIsims'
}
*************************************************************

local st = "$S_TIME"
local oper = "SAR Model for `MCsims' simulations"

*** Load the n-order contiguity matrices
spatwmat using "Data\W\W_`obs'rs.dta", name(W)
spatwmat using "Data\W\W2_`obs'rs.dta", name(W2_c)
spatwmat using "Data\W\W3_`obs'rs.dta", name(W3_c)	

*** Generate second- and third-order W matrix
mat W2 = W*W
mat W3 = W*W*W

tempvar id
gen `id' = _n

preserve
	use "Data\W\W_`obs'rs.dta", clear
	gen id = _n
	spmat dta Wspat V*, id(id) replace
restore

mat cons = J(`obs',1,1)

*** Provide access to the matrices in mata
mata: I = st_matrix("I")
mata: W = st_matrix("W")
mata: W2 = st_matrix("W2")
mata: W3 = st_matrix("W3")
mata: W2_c = st_matrix("W2_c")
mata: W3_c = st_matrix("W3_c")
mata: beta = st_numscalar("b1")
local b1 = b1

local stfull = "$S_TIME"
local operfull = "SAR Model for `MCsims' simulations"

*************************************************************
*** Calculate the True Effects
*************************************************************
tempname sartrue
postfile `sartrue' rho b1  /* 
*/	true_teffect true_deffect true_ieffect true_feffect true_heffect true_deffect_o0 true_teffect_o1 true_deffect_o1 true_ieffect_o1 true_teffect_o2 true_deffect_o2 true_ieffect_o2 true_teffect_o3 true_deffect_o3 true_ieffect_o3 /*
*/	using "Generate\SAR\SAR--True Effects.dta", replace every(20)

qui foreach r of numlist -0.8(.2)0.8 { 
	scalar rho = `r'
	mata: rho = st_numscalar("rho")

	matrix true = [`r', `b1']
	mata: true = st_matrix("true")
	mata: calceffsar(true, I, W, 1)
	
	scalar true_teffect = r(teffect)
	scalar true_deffect = r(deffect)
	scalar true_ieffect = r(ieffect)
	scalar true_feffect = r(feffect)
	scalar true_heffect = r(heffect)
	scalar true_deffect_o0 = r(deffect0)

	foreach o of numlist 1(1)3 {
		foreach e in t d i {
			scalar true_`e'effect_o`o' = r(`e'effect`o')
		}
	}	
		
	post `sartrue' (`r') (`b1') /* 
*/	(true_teffect) (true_deffect) (true_ieffect) (true_feffect) (true_heffect) (true_deffect_o0) (true_teffect_o1) (true_deffect_o1) (true_ieffect_o1) (true_teffect_o2) (true_deffect_o2) (true_ieffect_o2) (true_teffect_o3) (true_deffect_o3) (true_ieffect_o3)

}
	
postclose `sartrue'	

preserve
	use "Generate\SAR\SAR--True Effects.dta", clear
	sort rho
	save "Generate\SAR\SAR--True Effects.dta", replace
restore

*************************************************************
*** Monte Carlo Experiments
*************************************************************
tempname sar
postfile `sar' sim rho b1 str30 time b_x_sar se_x_sar rho_sar se_rho_sar double (rmse_sar rmse_xb_sar) /* 
*/	teffect teffect_se deffect deffect_se ieffect ieffect_se feffect feffect_se heffect heffect_se deffect_o0 deffect_se_o0 teffect_o1 teffect_se_o1 deffect_o1 deffect_se_o1 ieffect_o1 ieffect_se_o1 teffect_o2 teffect_se_o2 deffect_o2 deffect_se_o2 ieffect_o2 ieffect_se_o2 teffect_o3 teffect_se_o3 deffect_o3 deffect_se_o3 ieffect_o3 ieffect_se_o3 /*
*/	b_x_m1 se_x_m1 theta_1_m1 se_theta_1_m1 teffect_m1 teffect_se_m1 deffect_m1 deffect_se_m1 ieffect_m1 ieffect_se_m1 deffect_o0_m1 deffect_o0_se_m1 teffect_o1_m1 teffect_se_o1_m1 deffect_o1_m1 deffect_se_o1_m1 ieffect_o1_m1 ieffect_se_o1_m1 ar2_m1 vif_m1 /*
*/	b_x_m2 se_x_m2 theta_1_m2 se_theta_1_m2 theta_2_m2 se_theta_2_m2 teffect_m2 teffect_se_m2 deffect_m2 deffect_se_m2 ieffect_m2 ieffect_se_m2 feffect_m2 feffect_se_m2 heffect_m2 heffect_se_m2 deffect_o0_m2 deffect_se_o0_m2 teffect_o1_m2 teffect_se_o1_m2 deffect_o1_m2 deffect_se_o1_m2 ieffect_o1_m2 ieffect_se_o1_m2 teffect_o2_m2 teffect_se_o2_m2 deffect_o2_m2 deffect_se_o2_m2 ieffect_o2_m2 ieffect_se_o2_m2 p_spat_m2 ar2_m2 vif_m2 /*
*/	b_x_m3 se_x_m3 theta_1_m3 se_theta_1_m3 theta_2_m3 se_theta_2_m3 teffect_m3 teffect_se_m3 deffect_m3 deffect_se_m3 ieffect_m3 ieffect_se_m3 feffect_m3 feffect_se_m3 heffect_m3 heffect_se_m3 deffect_o0_m3 deffect_se_o0_m3 teffect_o1_m3 teffect_se_o1_m3 deffect_o1_m3 deffect_se_o1_m3 ieffect_o1_m3 ieffect_se_o1_m3 teffect_o2_m3 teffect_se_o2_m3 deffect_o2_m3 deffect_se_o2_m3 ieffect_o2_m3 ieffect_se_o2_m3 p_spat_m3 ar2_m3 vif_m3 /*
*/	b_x_m4 se_x_m4 theta_1_m4 se_theta_1_m4 theta_2_m4 se_theta_2_m4 theta_3_m4 se_theta_3_m4 teffect_m4 teffect_se_m4 deffect_m4 deffect_se_m4 ieffect_m4 ieffect_se_m4 feffect_m4 feffect_se_m4 heffect_m4 heffect_se_m4 deffect_o0_m4 deffect_se_o0_m4 teffect_o1_m4 teffect_se_o1_m4 deffect_o1_m4 deffect_se_o1_m4 ieffect_o1_m4 ieffect_se_o1_m4 teffect_o2_m4 teffect_se_o2_m4 deffect_o2_m4 deffect_se_o2_m4 ieffect_o2_m4 ieffect_se_o2_m4 teffect_o3_m4 teffect_se_o3_m4 deffect_o3_m4 deffect_se_o3_m4 ieffect_o3_m4 ieffect_se_o3_m4 p_spat_m4 ar2_m4 vif_m4 /*
*/	b_x_m5 se_x_m5 theta_1_m5 se_theta_1_m5 theta_2_m5 se_theta_2_m5 theta_3_m5 se_theta_3_m5 teffect_m5 teffect_se_m5 deffect_m5 deffect_se_m5 ieffect_m5 ieffect_se_m5 feffect_m5 feffect_se_m5 heffect_m5 heffect_se_m5 deffect_o0_m5 deffect_se_o0_m5 teffect_o1_m5 teffect_se_o1_m5 deffect_o1_m5 deffect_se_o1_m5 ieffect_o1_m5 ieffect_se_o1_m5 teffect_o2_m5 teffect_se_o2_m5 deffect_o2_m5 deffect_se_o2_m5 ieffect_o2_m5 ieffect_se_o2_m5 teffect_o3_m5 teffect_se_o3_m5 deffect_o3_m5 deffect_se_o3_m5 ieffect_o3_m5 ieffect_se_o3_m5 p_spat_m5 p_high_m5 ar2_m5 vif_m5  /*
*/	using "Generate\SAR\SAR.dta", replace every(20)

qui foreach r of numlist -0.8(.2)0.8 { 
	scalar rho = `r'
	mata: rho = st_numscalar("rho")
	
	qui foreach i of numlist 1(1)`MCsims' {
		local st1 = "$S_TIME"
		local oper1 = "Simulation `i'"
		
		*** Generate the error term
		drawnorm e, n(`obs') sd(`s')
		mkmat e, matrix(e) 

		*** Generate the X matrix
		generate x1 = -10+(10 - -10)*uniform()
		mkmat x1, matrix(X)

		*** Generate the spatial lags
		mat slag1 = W * X
		svmat slag1 

		mat slag2 = W2 * X 
		svmat slag2

		mat slag3 = W3 * X
		svmat slag3

		mat slag2_c = W2_c * X 
		svmat slag2_c

		mat slag3_c = W3_c * X
		svmat slag3_c
		
		*=========================================================
		*=========================================================
		*** Generate the SAR data-generating process (mata)

		mata: X = st_matrix("X")
		mata: e = st_matrix("e")
		mata: y = luinv(I-rho*W)*(X*beta) + luinv(I-rho*W)*(e) //generate the y from the multiplier and XB matrix

		getmata y // extract the y from mata to view

		*============================================================
		*============================================================	
		if `obs' > `QIsims' {
			set obs `obs'
}
		else {
			set obs `QIsims'
		}		

		*** Estimate the SAR
		spreg ml y x1 in 1/`obs', id(`id') dlmat(Wspat, eig)
		mat V = e(V)
		
		scalar b_x_sar = _b[y:x1]
		scalar se_x_sar = _se[y:x1]

		scalar rho_sar = _b[lambda:_cons]
		scalar se_rho_sar = _se[lambda:_cons]

		spreg_fit
		scalar rmse_sar = r(rmse)
		scalar rmse_xb_sar = r(rmse_xb)

		matrix coeff = [b_rho_sar, b_x_sar]
		matrix VCVM = [V[3,3], V[1,3] \ V[3,1], V[1,1]]
		
		preserve
			tempvar rho_draw beta_draw
			corr2data `rho_draw' `beta_draw', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `rho_draw' `beta_draw', matrix(draws)
		restore		
		
		mata: draws = st_matrix("draws")
		mata: calceffsar(draws, I, W, `QIsims')
		
		scalar teffect = r(teffect)
		scalar teffect_se = r(teffect_se)
		scalar deffect = r(deffect)
		scalar deffect_se = r(deffect_se)
		scalar ieffect = r(ieffect)
		scalar ieffect_se = r(ieffect_se)
		
		scalar feffect = r(feffect)
		scalar feffect_se = r(feffect_se)
		scalar heffect = r(heffect)
		scalar heffect_se = r(heffect_se)

		scalar deffect_o0 = r(deffect0)
		scalar deffect_se_o0 = r(deffect0_se)
		
		foreach o of numlist 1(1)3 {
			foreach e in t d i {
				scalar `e'effect_o`o' = r(`e'effect`o')
				scalar `e'effect_se_o`o' = r(`e'effect`o'_se)
			}
		}
		
		
		* Model 1: W 
		qui reg y x1 slag11
		mat b = e(b)
		mat V = e(V)
		matrix coeff = b[1,1..2]
		matrix VCVM = V[1..2,1..2]
		
		preserve		
			tempvar beta_m1 theta1_m1
			corr2data `beta_m1' `theta1_m1', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `beta_m1' `theta1_m1', matrix(draws_m1)	
		restore		
		
		mata: draws_m1 = st_matrix("draws_m1")
		mata: calceffslx1(draws_m1, W, `QIsims')
	
		foreach e in t d i {	
			scalar `e'effect_m1 = r(`e'effect)
			scalar `e'effect_se_m1 = r(`e'effect_se)
			foreach o of numlist 1(1)1 {
				scalar `e'effect_o`o'_m1 = r(`e'effect_o`o')
				scalar `e'effect_se_o`o'_m1 = r(`e'effect_se_o`o')
			}			
		}
	
		scalar deffect_o0_m1 = r(deffect_o0)
		scalar deffect_o0_se_m1 = r(deffect_se_o0)
	
		scalar b_x_m1 = _b[x1]
		scalar se_x_m1 = _se[x1]
	
		scalar b_theta_1_m1 = _b[slag11]
		scalar se_theta_1_m1 = _se[slag11]
	
		scalar ar2_m1 = e(r2_a)
		vif
		scalar vif_m1 = r(vif_1)
		
		
		
		* Model 2: W and W^2
		reg y x1 slag11 slag21
		mat b = e(b)
		mat V = e(V)
		matrix coeff = b[1,1..3]
		matrix VCVM = V[1..3,1..3]
		
		preserve		
			tempvar beta_m2 theta1_m2 theta2_m2
			corr2data `beta_m2' `theta1_m2' `theta2_m2', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `beta_m2' `theta1_m2' `theta2_m2', matrix(draws_m2)	
		restore		
		
		mata: draws_m2 = st_matrix("draws_m2")
		mata: calceffslx2(draws_m2, W, W2, `QIsims')

		foreach e in t d i {	
			scalar `e'effect_m2 = r(`e'effect)
			scalar `e'effect_se_m2 = r(`e'effect_se)
			foreach o of numlist 1(1)2 {
				scalar `e'effect_o`o'_m2 = r(`e'effect_o`o')
				scalar `e'effect_se_o`o'_m2 = r(`e'effect_se_o`o')
			}			
		}
		
		scalar feffect_m2 = r(feffect)
		scalar feffect_se_m2 = r(feffect_se)
		scalar heffect_m2 = r(heffect)
		scalar heffect_se_m2 = r(heffect_se)

		scalar deffect_o0_m2 = r(deffect_o0)
		scalar deffect_se_o0_m2 = r(deffect_se_o0)
		
		scalar b_x_m2 = _b[x1]
		scalar se_x_m2 = _se[x1]
	
		scalar b_theta_1_m2 = _b[slag11]
		scalar se_theta_1_m2 = _se[slag11]
	
		scalar b_theta_2_m2 = _b[slag21]
		scalar se_theta_2_m2 = _se[slag21]
	
		test slag11 = slag21 = 0
		scalar p_spat_m2 = r(p)
	
		scalar ar2_m2 = e(r2_a)
		vif
		scalar vif_m2 = r(vif_1)

		
		
		* Model 3: W and W2_c
		qui reg y x1 slag11 slag2_c1
		mat b = e(b)
		mat V = e(V)
		matrix coeff = b[1,1..3]
		matrix VCVM = V[1..3,1..3]
		
		preserve		
			tempvar beta_m3 theta1_m3 theta2_m3 
			corr2data `beta_m3' `theta1_m3' `theta2_m3', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `beta_m3' `theta1_m3' `theta2_m3', matrix(draws_m3)	
		restore		
		
		mata: draws_m3 = st_matrix("draws_m3")
		mata: calceffslx2(draws_m3, W, W2_c, `QIsims')

		foreach e in t d i {	
			scalar `e'effect_m3 = r(`e'effect)
			scalar `e'effect_se_m3 = r(`e'effect_se)
			foreach o of numlist 1(1)2 {
				scalar `e'effect_o`o'_m3 = r(`e'effect_o`o')
				scalar `e'effect_se_o`o'_m3 = r(`e'effect_se_o`o')
			}			
		}

		scalar feffect_m3 = r(feffect)
		scalar feffect_se_m3 = r(feffect_se)
		scalar heffect_m3 = r(heffect)
		scalar heffect_se_m3 = r(heffect_se)

		scalar deffect_o0_m3 = r(deffect_o0)
		scalar deffect_se_o0_m3 = r(deffect_se_o0)
		
		scalar b_x_m3 = _b[x1]
		scalar se_x_m3 = _se[x1]
	
		scalar b_theta_1_m3 = _b[slag11]
		scalar se_theta_1_m3 = _se[slag11]
	
		scalar b_theta_2_m3 = _b[slag2_c1]
		scalar se_theta_2_m3 = _se[slag2_c1]
	
		test slag11 = slag2_c1 = 0
		scalar p_spat_m3 = r(p)
	
		scalar ar2_m3 = e(r2_a)
		vif
		scalar vif_m3 = r(vif_1)

		
		
		* Model 4: W, W^2 and W^3
		reg y x1 slag11 slag21 slag31
		mat b = e(b)
		mat V = e(V)
		matrix coeff = b[1,1..4]
		matrix VCVM = V[1..4,1..4]
		
		preserve		
			tempvar beta_m4 theta1_m4 theta2_m4 theta3_m4
			corr2data `beta_m4' `theta1_m4' `theta2_m4' `theta3_m4', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `beta_m4' `theta1_m4' `theta2_m4' `theta3_m4', matrix(draws_m4)	
		restore
		
		mata: draws_m4 = st_matrix("draws_m4")
		mata: calceffslx3(draws_m4, W, W2, W3, `QIsims')		
		
		foreach e in t d i {
			scalar `e'effect_m4 = r(`e'effect)
			scalar `e'effect_se_m4 = r(`e'effect_se)
						
			foreach o of numlist 1(1)3 {
				scalar `e'effect_o`o'_m4 = r(`e'effect_o`o')
				scalar `e'effect_se_o`o'_m4 = r(`e'effect_se_o`o')
			}
		}		
		
		scalar feffect_m4 = r(feffect)
		scalar feffect_se_m4 = r(feffect_se)
		scalar heffect_m4 = r(heffect)
		scalar heffect_se_m4 = r(heffect_se)
		
		scalar deffect_o0_m4 = r(deffect_o0)
		scalar deffect_se_o0_m4 = r(deffect_se_o0)
		
		scalar b_x_m4 = _b[x1]
		scalar se_x_m4 = _se[x1]
	
		scalar b_theta_1_m4 = _b[slag11]
		scalar se_theta_1_m4 = _se[slag11]
	
		scalar b_theta_2_m4 = _b[slag21]
		scalar se_theta_2_m4 = _se[slag21]

		scalar b_theta_3_m4 = _b[slag31]
		scalar se_theta_3_m4 = _se[slag31]
		
		test slag11 = slag21 = slag31 = 0
		scalar p_spat_m4 = r(p)
	
		scalar ar2_m4 = e(r2_a)
		vif
		scalar vif_m4 = r(vif_1)
	
	
	
		* Model 5: W, W2_c and W3_c
		qui reg y x1 slag11 slag2_c1 slag3_c1
		mat b = e(b)
		mat V = e(V)
		matrix coeff = b[1,1..4]
		matrix VCVM = V[1..4,1..4]
		
		preserve		
			tempvar beta_m5 theta1_m5 theta2_m5 theta3_m5
			corr2data `beta_m5' `theta1_m5' `theta2_m5' `theta3_m5', n(`QIsims') means(coeff) cov(VCVM) clear
			mkmat `beta_m5' `theta1_m5' `theta2_m5' `theta3_m5', matrix(draws_m5)
		restore		
		
		mata: draws_m5 = st_matrix("draws_m5")
		mata: calceffslx3(draws_m5, W, W2_c, W3_c, `QIsims')

		foreach e in t d i {
			scalar `e'effect_m5 = r(`e'effect)
			scalar `e'effect_se_m5 = r(`e'effect_se)
						
			foreach o of numlist 1(1)3 {
				scalar `e'effect_o`o'_m5 = r(`e'effect_o`o')
				scalar `e'effect_se_o`o'_m5 = r(`e'effect_se_o`o')
			}
		}		
		
		scalar feffect_m5 = r(feffect)
		scalar feffect_se_m5 = r(feffect_se)
		scalar heffect_m5 = r(heffect)
		scalar heffect_se_m5 = r(heffect_se)
		
		scalar deffect_o0_m5 = r(deffect_o0)
		scalar deffect_se_o0_m5 = r(deffect_se_o0)
		
		scalar b_x_m5 = _b[x1]
		scalar se_x_m5 = _se[x1]
	
		scalar b_theta_1_m5 = _b[slag11]
		scalar se_theta_1_m5 = _se[slag11]
	
		scalar b_theta_2_m5 = _b[slag2_c1]
		scalar se_theta_2_m5 = _se[slag2_c1]
	
		scalar b_theta_3_m5 = _b[slag3_c1]
		scalar se_theta_3_m5 = _se[slag3_c1]
	
		test slag11 = slag2_c1 = slag3_c1 = 0
		scalar p_spat_m5 = r(p)
	
		test slag2_c1 = slag3_c1 = 0
		scalar p_high_m5 = r(p)
	
		scalar ar2_m5 = e(r2_a)
		vif
		scalar vif_m5 = r(vif_1)
		
		if mod(`i',5) == 0 {
			nois display "." _c
			if mod(`i',100) == 0 {
				nois display "$S_DATE : `c(current_time)': Simulation #`i'"
			}
		}		

		elapse "`st1'" "`oper1'"		
		
		post `sar' (`i') (`r') (`b1') ("$S_elap") (b_x_sar) (se_x_sar) (rho_sar) (se_rho_sar) (rmse_sar) (rmse_xb_sar) /* 
		*/	(teffect) (teffect_se) (deffect) (deffect_se) (ieffect) (ieffect_se) (feffect) (feffect_se) (heffect) (heffect_se) (deffect_o0) (deffect_se_o0) (teffect_o1) (teffect_se_o1) (deffect_o1) (deffect_se_o1) (ieffect_o1) (ieffect_se_o1) (teffect_o2) (teffect_se_o2) (deffect_o2) (deffect_se_o2) (ieffect_o2) (ieffect_se_o2) (teffect_o3) (teffect_se_o3) (deffect_o3) (deffect_se_o3) (ieffect_o3) (ieffect_se_o3) /*
		*/	(b_x_m1) (se_x_m1) (b_theta_1_m1) (se_theta_1_m1) (teffect_m1) (teffect_se_m1) (deffect_m1) (deffect_se_m1) (ieffect_m1) (ieffect_se_m1) (deffect_o0_m1) (deffect_o0_se_m1) (teffect_o1_m1) (teffect_se_o1_m1) (deffect_o1_m1) (deffect_se_o1_m1) (ieffect_o1_m1) (ieffect_se_o1_m1) (ar2_m1) (vif_m1) /*
		*/	(b_x_m2) (se_x_m2) (b_theta_1_m2) (se_theta_1_m2) (b_theta_2_m2) (se_theta_2_m2) (teffect_m2) (teffect_se_m2) (deffect_m2) (deffect_se_m2) (ieffect_m2) (ieffect_se_m2) (feffect_m2) (feffect_se_m2) (heffect_m2) (heffect_se_m2) (deffect_o0_m2) (deffect_se_o0_m2) (teffect_o1_m2) (teffect_se_o1_m2) (deffect_o1_m2) (deffect_se_o1_m2) (ieffect_o1_m2) (ieffect_se_o1_m2) (teffect_o2_m2) (teffect_se_o2_m2) (deffect_o2_m2) (deffect_se_o2_m2) (ieffect_o2_m2) (ieffect_se_o2_m2) (p_spat_m2) (ar2_m2) (vif_m2) /*
		*/	(b_x_m3) (se_x_m3) (b_theta_1_m3) (se_theta_1_m3) (b_theta_2_m3) (se_theta_2_m3) (teffect_m3) (teffect_se_m3) (deffect_m3) (deffect_se_m3) (ieffect_m3) (ieffect_se_m3) (feffect_m3) (feffect_se_m3) (heffect_m3) (heffect_se_m3) (deffect_o0_m3) (deffect_se_o0_m3) (teffect_o1_m3) (teffect_se_o1_m3) (deffect_o1_m3) (deffect_se_o1_m3) (ieffect_o1_m3) (ieffect_se_o1_m3) (teffect_o2_m3) (teffect_se_o2_m3) (deffect_o2_m3) (deffect_se_o2_m3) (ieffect_o2_m3) (ieffect_se_o2_m3) (p_spat_m3) (ar2_m3) (vif_m3) /*
		*/	(b_x_m4) (se_x_m4) (b_theta_1_m4) (se_theta_1_m4) (b_theta_2_m4) (se_theta_2_m4) (b_theta_3_m4) (se_theta_3_m4) (teffect_m4) (teffect_se_m4) (deffect_m4) (deffect_se_m4) (ieffect_m4) (ieffect_se_m4) (feffect_m4) (feffect_se_m4) (heffect_m4) (heffect_se_m4) (deffect_o0_m4) (deffect_se_o0_m4) (teffect_o1_m4) (teffect_se_o1_m4) (deffect_o1_m4) (deffect_se_o1_m4) (ieffect_o1_m4) (ieffect_se_o1_m4) (teffect_o2_m4) (teffect_se_o2_m4) (deffect_o2_m4) (deffect_se_o2_m4) (ieffect_o2_m4) (ieffect_se_o2_m4) (teffect_o3_m4) (teffect_se_o3_m4) (deffect_o3_m4) (deffect_se_o3_m4) (ieffect_o3_m4) (ieffect_se_o3_m4) (p_spat_m4) (ar2_m4) (vif_m4) /*
		*/	(b_x_m5) (se_x_m5) (b_theta_1_m5) (se_theta_1_m5) (b_theta_2_m5) (se_theta_2_m5) (b_theta_3_m5) (se_theta_3_m5) (teffect_m5) (teffect_se_m5) (deffect_m5) (deffect_se_m5) (ieffect_m5) (ieffect_se_m5) (feffect_m5) (feffect_se_m5) (heffect_m5) (heffect_se_m5) (deffect_o0_m5) (deffect_se_o0_m5) (teffect_o1_m5) (teffect_se_o1_m5) (deffect_o1_m5) (deffect_se_o1_m5) (ieffect_o1_m5) (ieffect_se_o1_m5) (teffect_o2_m5) (teffect_se_o2_m5) (deffect_o2_m5) (deffect_se_o2_m5) (ieffect_o2_m5) (ieffect_se_o2_m5) (teffect_o3_m5) (teffect_se_o3_m5) (deffect_o3_m5) (deffect_se_o3_m5) (ieffect_o3_m5) (ieffect_se_o3_m5) (p_spat_m5) (p_high_m5) (ar2_m5) (vif_m5)

		keep in 1/`obs'	
		drop e x1 y slag1* slag2* slag3*
		mata: mata drop draws draws_m1 draws_m2 draws_m3 draws_m4 draws_m5 /* fbeffect_m3 fbeffect_m5	*/
	}
}
	
postclose `sar'	
elapse "`stfull'" "`operfull'"
